home *** CD-ROM | disk | FTP | other *** search
/ Delphi 5 for Professionals / DELPHI5.iso / AddOns / Components / TEECHART / Src Code / TEEPOLY.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1998-10-24  |  4.7 KB  |  179 lines
  1. {****************************************************}
  2. {   TeeChart Pro 4.0                     }
  3. {   Polynomic routines                               }
  4. {   Copyright (c) 1995-1998 by David Berneda         }
  5. {****************************************************}
  6. Unit TeePoly;
  7.  
  8. interface
  9.  
  10. Uses Classes;
  11.  
  12. Const MaxPolyDegree=20;    { maximum number of polynomial degree }
  13.       MaxPolyPoints=16300; { maximum number of points }
  14.  
  15. Type Float         = Extended;
  16.      PFloat        = ^Float;
  17.      TDegreeVector = Array[1..MaxPolyDegree] of Float;
  18.      TMatrix       = Array[1..MaxPolyDegree,1..MaxPolyDegree] of Float;
  19.  
  20.      PVector       = ^TVector;
  21.      TVector       = Array[0..MaxPolyPoints] of PFloat;
  22.  
  23. Function CalcFitting(PolyDegree:Integer; Const Answer:TDegreeVector; Const XWert:Float):Float;
  24. Procedure PolyFitting( NumPoints:Longint; PolyDegree:Integer;
  25.                        X,Y:PVector; var Answer:TDegreeVector);
  26.  
  27. Procedure SetVectorValue(Const V:TVector; Index:Longint; Const Value:Float);
  28. Function GetVectorValue(Const V:TVector; Index:Longint):Float;
  29.  
  30. implementation
  31.  
  32. Uses SysUtils,TeeProCo;
  33.  
  34. Function GetVectorValue(Const V:TVector; Index:Longint):Float;
  35. Begin
  36.   result:=V[Index]^;
  37. end;
  38.  
  39. Procedure SetVectorValue(Const V:TVector; Index:Longint; Const Value:Float);
  40. begin
  41.   V[Index]^:=Value;
  42. end;
  43.  
  44. Function GaussianFitting( NumDegree:Integer;
  45.                           Var M:TMatrix;
  46.                           Var Y,X:TDegreeVector;
  47.                           Const Error:Float):Float;
  48. var i        : Integer;
  49.     j        : Integer;
  50.     k        : Integer;
  51.     MaxIndex : Integer;
  52.     Change   : Integer;
  53.     MaxEl    : Float;
  54.     Temp     : Float;
  55.     Wert     : Float;
  56. begin
  57.   Change:=0;
  58.   for i:=1 to NumDegree-1 do
  59.   begin
  60.     MaxIndex:=i;
  61.     MaxEl:=Abs(M[i,i]);
  62.     for j:=i+1 to NumDegree do
  63.     if Abs(M[j,i]) > MaxEl then
  64.     begin
  65.       MaxEl:=Abs(M[j,i]);
  66.       MaxIndex:=j;
  67.     end;
  68.     if MaxIndex <> i then
  69.     begin
  70.        for j:=i to NumDegree do
  71.        begin
  72.          Temp:=M[i,j];
  73.          M[i,j]:=M[MaxIndex,j];
  74.          M[MaxIndex,j]:=Temp;
  75.        end;
  76.        Temp:=Y[i];
  77.        Y[i]:=Y[MaxIndex];
  78.        Y[MaxIndex]:=Temp;
  79.        Inc(Change);
  80.     end;
  81.     if Abs(M[i,i]) < Error then
  82.     Begin
  83.       result:=0;
  84.       exit;
  85.     end;
  86.     for j:=i+1 to NumDegree do
  87.     begin
  88.       Wert:=M[j,i]/M[i,i];
  89.       for k:=i+1 to NumDegree do M[j,k]:=M[j,k]-Wert*M[i,k];
  90.       Y[j]:=Y[j]-Wert*Y[i];
  91.     end;
  92.   end;
  93.   if Abs(M[NumDegree,NumDegree])< Error then
  94.   Begin
  95.     Result:=0;
  96.     Exit;
  97.   end;
  98.   Result:=1;
  99.   for i:=NumDegree DownTo 1 do
  100.   begin
  101.     Wert:=0;
  102.     for j:=i+1 to NumDegree do Wert:=Wert + M[i,j]*X[j];
  103.     X[i]:=(Y[i]-Wert)/M[i,i];
  104.     result:=result*M[i,i];
  105.   end;
  106.   if Odd(Change) then result:=-result;
  107. end;
  108.  
  109. Procedure FKT(PolyDegree:Integer; Const xarg:Float; Var PHI:TDegreeVector);
  110. var t:Integer;
  111. begin
  112.   PHI[1]:=1;
  113.   for t:=2 to PolyDegree do PHI[t]:=PHI[t-1]*xarg;
  114. end;
  115.  
  116. Function CalcFitting( PolyDegree:Integer;
  117.                       Const Answer:TDegreeVector;
  118.                       Const xwert:Float):Float;
  119. var PHI : TDegreeVector;
  120.     t   : Integer;
  121. begin
  122.   result:=0;
  123.   FKT(PolyDegree,xwert,PHI);
  124.   for t:=1 to PolyDegree do result:=result+Answer[t]*PHI[t];
  125. end;
  126.  
  127. Procedure PolyFitting( NumPoints:Longint; PolyDegree:Integer; X,Y:PVector;
  128.                        var Answer:TDegreeVector);
  129. var  t     : Integer;
  130.      tt    : Integer;
  131.      l     : Integer;
  132.      PHI   : TDegreeVector;
  133.      B     : TDegreeVector;
  134.      F     : Array[1..MaxPolyDegree] of PVector;
  135.      M     : TMatrix;
  136. begin
  137.   for t:=1 to PolyDegree do F[t]:=nil;
  138.   for t:=1 to PolyDegree do
  139.   Begin
  140.     New(F[t]);
  141.     for tt:=0 to NumPoints do F[t]^[tt]:=New(PFloat);
  142.   end;
  143.   try
  144.     for t:=1 to NumPoints do         { Prepare the approximation }
  145.     begin
  146.       FKT(PolyDegree,GetVectorValue(X^,t),PHI);
  147.       for tt:=1 to PolyDegree do SetVectorValue(F[tt]^,t,PHI[tt]);
  148.     end;
  149.     for tt:=1 to PolyDegree  do         { Build the matrix of the LinEqu. }
  150.     for t:=1 to PolyDegree  do
  151.     begin
  152.       M[t,tt]:=0;
  153.       for l:=1 to NumPoints do
  154.           M[t,tt]:=M[t,tt]+GetVectorValue(F[t]^,l)*GetVectorValue(F[tt]^,l);
  155.       M[tt,t]:=M[t,tt];
  156.     end;
  157.     for t:=1 to PolyDegree  do
  158.     begin
  159.       B[t]:=0;
  160.       for l:=1 to NumPoints do
  161.           B[t]:=B[t]+GetVectorValue(F[t]^,l)*GetVectorValue(Y^,l);
  162.     end;
  163.  
  164.     if GaussianFitting(PolyDegree,M,B,Answer,1.0e-15)=0 then
  165.        Raise Exception.Create(TeeMsg_FittingError);
  166.   finally
  167.     for t:=1 to PolyDegree do
  168.     if Assigned(F[t]) then
  169.     Begin
  170.       for tt:=0 to NumPoints do Dispose(PFloat(F[t]^[tt]));
  171.       Dispose(F[t]);
  172.     end;
  173.   end;
  174. end;
  175.  
  176. end.
  177.  
  178.  
  179.